# Load necessary librarieslibrary(forecast)library(tseries)library(tidyverse)library(tidymodels)library(kableExtra)library(modeltime)library(timetk)library(cowplot)# setting the seedset.seed(123)# Loading Data ----------------------------------------------------------------daily_deaths <-read_csv("data/new_daily_deaths.csv") %>%select(daily_deaths)east <-read_csv("data/east_daily.csv")midwest <-read_csv("data/midwest_daily.csv")south <-read_csv("data/south_daily.csv")west <-read_csv("data/west_daily.csv")
2. Data Splitting
Code
# Data Splitting ---------------------------------------------------------------east_splits <-time_series_split(east, initial ="1022 days", assess ="114 days")east_train <-training(east_splits)east_test <-testing(east_splits)midwest_splits <-time_series_split(midwest, initial ="1022 days", assess ="114 days")midwest_train <-training(midwest_splits)midwest_test <-testing(midwest_splits)south_splits <-time_series_split(south, initial ="1022 days", assess ="114 days")south_train <-training(south_splits)south_test <-testing(south_splits)west_splits <-time_series_split(west, initial ="1022 days", assess ="114 days")west_train <-training(west_splits)west_test <-testing(west_splits)## Convert Into Time Series ----------------------------------------------------east_train_ts <-ts(east_train %>%select(daily_deaths))east_test_ts <-ts(east_test %>%select(daily_deaths))east_ts <-ts(east %>%select(daily_deaths))midwest_train_ts <-ts(midwest_train %>%select(daily_deaths))midwest_test_ts <-ts(midwest_test %>%select(daily_deaths))midwest_ts <-ts(midwest %>%select(daily_deaths))south_train_ts <-ts(south_train %>%select(daily_deaths))south_test_ts <-ts(south_test %>%select(daily_deaths))south_ts <-ts(south %>%select(daily_deaths))west_train_ts <-ts(west_train %>%select(daily_deaths))west_test_ts <-ts(west_test %>%select(daily_deaths))west_ts <-ts(west %>%select(daily_deaths))## Plotting Training and Testing -----------------------------------------------par(mfrow=c(2,2))ggplot() +geom_line(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = east_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = east_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="East Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="Midwest Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = south_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = south_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="South Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = west_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = west_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="West Region Training vs. Testing Data")
3. Testing for Stationarity
Code
### Using Augmented Dickey-Fuller Test -----------------------------------------adf_east <-adf.test(east_ts)print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary
Augmented Dickey-Fuller Test
data: east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_midwest <-adf.test(midwest_ts)print(adf_midwest) # p-value is 0.325 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary
Code
midwest_ts <-diff(midwest_ts)adf_midwest <-adf.test(midwest_ts)print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_south <-adf.test(south_ts)print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary
Code
south_ts <-diff(south_ts)adf_south <-adf.test(south_ts)print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_west <-adf.test(west_ts)print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary
Code
west_ts <-diff(west_ts)adf_west <-adf.test(west_ts)print(adf_west) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
## ACF & PACF Plots ------------------------------------------------------------par(mfrow=c(3, 1))acf(east_ts, main ="East")pacf(east_ts, main ="East")plot(east_ts, main ="Daily Deaths: East")
Code
acf(midwest_ts, main ="Midwest")pacf(midwest_ts, main ="Midwest")plot(midwest_ts, main ="Daily Deaths: Midwest")
Code
acf(south_ts, main ="South")pacf(south_ts, main ="South")plot(south_ts, main ="Daily Deaths: South")
Code
acf(west_ts, main ="West")pacf(west_ts, main ="West")plot(south_ts, main ="Daily Deaths: South")
Code
plot(west_ts, main ="Daily Deaths: West")
5. Time Series Decomposition
Code
# total daily death decompositiondaily_ts <-ts(daily_deaths, frequency =365)decomposed_ts <-decompose(daily_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(daily_ts, main ="Total Original Time Series", ylab ="Value", col ="blue")plot(decomposed_ts$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_ts$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_ts$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
# East Decomposition -----------------------------------------------------------east_daily <- east %>%select(daily_deaths)east_ts <-ts(east_daily, frequency =365)decomposed_east <-decompose(east_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(east_ts, main ="East Original Time Series", ylab ="Value", col ="blue")plot(decomposed_east$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_east$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_east$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## Midwest Decomposition ----------------------------------------------------------midwest_daily <- midwest %>%select(daily_deaths)midwest_ts <-ts(midwest_daily, frequency =365)decomposed_midwest <-decompose(midwest_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(midwest_ts, main ="Midwest Original Time Series", ylab ="Value", col ="blue")plot(decomposed_midwest$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_midwest$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_midwest$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## South Decomposition ----------------------------------------------------------south_daily <- south %>%select(daily_deaths)south_ts <-ts(south_daily, frequency =365)decomposed_south <-decompose(south_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(south_ts, main ="South Original Time Series", ylab ="Value", col ="blue")plot(decomposed_south$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_south$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_south$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## West Decomposition ----------------------------------------------------------west_daily <- west %>%select(daily_deaths)west_ts <-ts(west_daily, frequency =365)decomposed_west <-decompose(west_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(west_ts, main ="West Original Time Series", ylab ="Value", col ="blue")plot(decomposed_west$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_west$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_west$random, main ="Residual Component", ylab ="Value", col ="purple")
6. ARIMA Model
6.1 Initial Model
Code
# Building ARIMA Model ---------------------------------------------------------east_arima <-arima(east_train_ts, order=c(1,0,1))midwest_arima <-arima(midwest_train_ts, order=c(1,1,1))south_arima <-arima(south_train_ts, order=c(1,1,1))west_arima <-arima(west_train_ts, order=c(1,1,1))## Summary & Forecast ---------------------------------------------------------summary(east_arima)
Call:
arima(x = east_train_ts, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
0.9176 -0.0388 197.3738
s.e. 0.0147 0.0493 39.7998
sigma^2 estimated as 12153: log likelihood = -6257.17, aic = 12522.33
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1151898 110.2407 64.66739 -Inf Inf 0.9918341 0.008515868
Code
summary(midwest_arima)
Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.2237 -0.9343
s.e. 0.0313 0.0090
sigma^2 estimated as 16650: log likelihood = -6405.86, aic = 12817.73
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7797006 128.9735 78.3255 NaN Inf 0.7410171 -0.06964661
Code
summary(south_arima)
Call:
arima(x = south_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.1913 -0.9453
s.e. 0.0317 0.0095
sigma^2 estimated as 44608: log likelihood = -6908.51, aic = 13823.02
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.277928 211.1023 129.556 NaN Inf 0.7498262 -0.01884038
Code
summary(west_arima)
Call:
arima(x = west_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.0842 -0.9603
s.e. 0.0319 0.0072
sigma^2 estimated as 10324: log likelihood = -6162.22, aic = 12330.44
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6807731 101.558 62.19197 -Inf Inf 0.7886397 -0.009430079
Code
east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))
### South Regionsouth_initial_arima_acc <- forecast::accuracy(south_values)south_initial_arima_acc %>%kbl(caption ="<center>Initial South Accuracy<center>") %>%kable_styling()
### West Regionwest_initial_arima_acc <- forecast::accuracy(west_values)west_initial_arima_acc %>%kbl(caption ="<center>Initial West Accuracy<center>") %>%kable_styling()
# Grid Search for p and q ------------------------------------------------------p_values <-0:5# AR orderd_values <-0:0# I orderq_values <-0:5# MA order## East - Perform grid search for ARIMA parameters ----------------------------best_east_model <-NULLbest_east_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(east_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_east_aic) { best_east_model <- arima_model best_east_aic <-AIC(arima_model) best_east_order <-c(p, d, q) } } } }}
Code
## Midwest - Perform grid search for ARIMA parameters -------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_midwest_model <-NULLbest_midwest_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(midwest_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_midwest_aic) { best_midwest_model <- arima_model best_midwest_aic <-AIC(arima_model) best_midwest_order <-c(p, d, q) } } } }}
Code
## South - Perform grid search for ARIMA parameters ---------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_south_model <-NULLbest_south_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(south_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_south_aic) { best_south_model <- arima_model best_south_aic <-AIC(arima_model) best_south_order <-c(p, d, q) } } } }}
Code
## West - Perform grid search for ARIMA parameters ----------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_west_model <-NULLbest_west_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(west_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_west_aic) { best_west_model <- arima_model best_west_aic <-AIC(arima_model) best_west_order <-c(p, d, q) } } } }}
### South Regionsouth_tuned_arima_acc <- forecast::accuracy(south_values)south_tuned_arima_acc %>%kbl(caption ="<center>Tuned South Accuracy<center>") %>%kable_styling()
### West Regionwest_tuned_arima_acc <- forecast::accuracy(west_values)west_tuned_arima_acc %>%kbl(caption ="<center>Tuned West Accuracy<center>") %>%kable_styling()
# Building SARIMA model ------------------------------------------------------east_sarima <-arima(east_train_ts, seasonal =list(order =c(1,1,1)))midwest_sarima <-arima(midwest_train_ts, seasonal =list(order =c(1,1,1)))south_sarima <-arima(south_train_ts, seasonal =list(order =c(1,1,1)))west_sarima <-arima(west_train_ts, seasonal =list(order =c(1,1,1)))# Summary of the SARIMA modelsummary(east_sarima)
Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
0.4706 -0.8016
s.e. 0.0382 0.0211
sigma^2 estimated as 11156: log likelihood = -6206.64, aic = 12419.27
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1606608 105.5715 61.2767 -Inf Inf 0.9398295 0.1140208
Code
summary(midwest_sarima)
Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.2237 -0.9343
s.e. 0.0313 0.0090
sigma^2 estimated as 16650: log likelihood = -6405.86, aic = 12817.73
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7797006 128.9735 78.3255 NaN Inf 0.7410171 -0.06964661
Code
summary(south_sarima)
Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.1913 -0.9453
s.e. 0.0317 0.0095
sigma^2 estimated as 44608: log likelihood = -6908.51, aic = 13823.02
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.277928 211.1023 129.556 NaN Inf 0.7498262 -0.01884038
Code
summary(west_sarima)
Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.0842 -0.9603
s.e. 0.0319 0.0072
sigma^2 estimated as 10324: log likelihood = -6162.22, aic = 12330.44
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6807731 101.558 62.19197 -Inf Inf 0.7886397 -0.009430079
Code
## Forecasting with the SARIMA model -------------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
7.2 Accuracy
Code
## Plotting & Accuracy## East Regioneast_initial_sarima_acc <- forecast::accuracy(east_values)east_initial_sarima_acc %>%kbl(caption ="<center>Initial SARIMA East Accuracy<center>") %>%kable_styling()
## Define the parameter grids --------------------------------------------------p_grid <-0:5# AR parameterd_grid <-0:0# Differencing parameterq_grid <-0:5# MA parameterP_grid <-0:1# Seasonal AR parameterD_grid <-0:1# Seasonal differencing parameterQ_grid <-0:1# Seasonal MA parameters <-12# Seasonal period## East Grid Search ------------------------------------------------------------# Create East Regioneast_train <- east_train %>%select(daily_deaths)# Create an empty data frame to store the resultseast_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(east_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame east_results <-rbind(east_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}d_grid <-0:1# Differencing parameter## Midwest Grid Search ------------------------------------------------------------# Create Midwest Regionmidwest_train <- midwest_train %>%select(daily_deaths)# Create an empty data frame to store the resultsmidwest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(midwest_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame midwest_results <-rbind(midwest_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## South Grid Search ------------------------------------------------------------# Create South Regionsouth_train <- south_train %>%select(daily_deaths)# Create an empty data frame to store the resultssouth_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(south_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame south_results <-rbind(south_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## West Grid Search ------------------------------------------------------------# Create West Regionwest_train <- west_train %>%select(daily_deaths)# Create an empty data frame to store the resultswest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(west_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame west_results <-rbind(west_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}
# Rebuilding SARIMA Models -----------------------------------------------------east_sarima <-arima(east_train_ts, order =c(5, 0, 4), seasonal =list(order =c(0, 1, 1), period =12))midwest_sarima <-arima(midwest_train_ts, order =c(5, 1, 5),seasonal =list(order =c(0, 0, 1), period =12))south_sarima <-arima(south_train_ts, order =c(5, 1, 5),seasonal =list(order =c(0, 0, 0), period =12))west_sarima <-arima(west_train_ts, order =c(5, 0, 5),seasonal =list(order =c(1, 1, 0), period =12))## Forecasting with the NEW SARIMA model ---------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
Code
## Plotting & Accuracy## East Regioneast_tuned_sarima_acc <- forecast::accuracy(east_values)east_tuned_sarima_acc %>%kbl(caption ="<center>Tuned SARIMA East Accuracy<center>") %>%kable_styling()
# Building AutoARIMA Model -----------------------------------------------------east_arima <-auto.arima(east_train_ts, seasonal =TRUE)midwest_arima <-auto.arima(midwest_train_ts, seasonal =TRUE)south_arima <-auto.arima(south_train_ts, seasonal =TRUE)west_arima <-auto.arima(west_train_ts, seasonal =TRUE)# Summary of the ARIMA modelsummary(east_arima)
Series: east_train_ts
ARIMA(5,1,2)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2
0.5142 -0.7304 -0.1394 -0.3479 -0.2207 -1.0295 0.7353
s.e. 0.0480 0.0344 0.0426 0.0329 0.0418 0.0410 0.0265
sigma^2 = 6861: log likelihood = -5956.39
AIC=11928.79 AICc=11928.93 BIC=11968.22
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1506172 82.50546 47.1043 NaN Inf 0.7224607 0.01707355
Code
summary(midwest_arima)
Series: midwest_train_ts
ARIMA(5,1,4)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
-0.0349 -0.7605 -0.0404 -0.1320 -0.4194 -1.1696 0.9929 -0.9722
s.e. 0.0461 0.0370 0.0572 0.0371 0.0377 0.0441 0.0492 0.0459
ma4
0.4244
s.e. 0.0381
sigma^2 = 10525: log likelihood = -6168.66
AIC=12357.32 AICc=12357.54 BIC=12406.6
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.385353 102.0896 59.56282 NaN Inf 0.5635082 -0.05532994
Code
summary(south_arima)
Series: south_train_ts
ARIMA(1,1,2)
Coefficients:
ar1 ma1 ma2
0.4799 -1.6879 0.7236
s.e. 0.0578 0.0428 0.0405
sigma^2 = 42636: log likelihood = -6884.08
AIC=13776.16 AICc=13776.2 BIC=13795.87
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.7155848 206.0793 126.7407 NaN Inf 0.7335323 0.008694211
Code
summary(west_arima)
Series: west_train_ts
ARIMA(5,1,3)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
0.4572 -0.6143 -0.2778 -0.2654 -0.2584 -1.8532 1.6316 -0.6651
s.e. 0.0588 0.0355 0.0387 0.0340 0.0435 0.0544 0.0868 0.0458
sigma^2 = 6300: log likelihood = -5907.95
AIC=11833.91 AICc=11834.09 BIC=11878.26
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.44309 79.02328 45.78169 NaN Inf 0.5805454 0.003263827
Code
## Forecasting Test Data -------------------------------------------------------east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))
8.2 Accuracy
Code
## Plotting & Accuracy## East Regioneast_auto_acc <- forecast::accuracy(east_values)east_auto_acc %>%kbl(caption ="<center>Auto ARIMA East Accuracy<center>") %>%kable_styling()